#include "HY1C_out.h"
#include <sys/types.h>
#include <unistd.h>

/* ============================================================================ */
/* module sstref.c - retrieve sst reference temperature from ancillary file     */
/*                                                                              */
/* Written By: B. Franz, NASA/SIMBIOS, August 2003.                             */
/*                                                                              */
/* ============================================================================ */

#include "l12_proto.h"
#include "hdf_utils.h"
#include <netcdf.h>

static float sstbad = BAD_FLT;
static int32_t  format = -1;


/* ----------------------------------------------------------------------------------- */
/* get_avhrrsst() - read and interpolate Reynolds 0.25-deg netcdf OI files             */
/*                                                                                     */
/* B. Franz, SAIC, July 2008.                                                          */
/* ----------------------------------------------------------------------------------- */

#define OI4NX 1440
#define OI4NY 720

float get_avhrrsst(char *sstfile, float lon, float lat)
{
    static int   firstCall = 1;
    static int   nx = OI4NX;
    static int   ny = OI4NY;
    static float dx = 360.0/OI4NX;
    static float dy = 180.0/OI4NY;
    static float sstref[OI4NY+2][OI4NX+2];

    float sst = sstbad;
    int   i,j,ii;
    int   ntmp;
    float xx,yy;
    float t,u,w[4],wt;
    float reftmp[2][2];
    float ftmp;

    if (firstCall) {
	typedef int16 ssttmp_t[OI4NX];
	ssttmp_t *ssttmp;
	ssttmp = (ssttmp_t*) malloc(OI4NY*sizeof(ssttmp_t));

//        int16 ssttmp[OI4NY][OI4NX];

        char *tmp_str;
        char  name   [H4_MAX_NC_NAME]  = "";
        char  sdsname[H4_MAX_NC_NAME]  = "";
        int ncid, grpid, ndims, nvars, ngatts, unlimdimid;
        int32 sd_id;
        int32 sds_id; 
        int32 rank; 
        int32 nt; 
        int32 dims[H4_MAX_VAR_DIMS]; 
        int32 nattrs;
        int32 start[4] = {0,0,0,0}; 
        int32 status;
        float slope;
        float offset;

        firstCall = 0;

        HY1C_out("Loading 0.25-deg OI Reynolds SST reference from %s\n",sstfile);
        HY1C_out("\n");

        /* Open the file */
	strcpy(sdsname,"sst");

	/* try netCDF first */
        if (nc_open(sstfile, NC_NOWRITE, &ncid) == 0) { 
		
            status = nc_inq(ncid, &ndims, &nvars, &ngatts, &unlimdimid);
            status = nc_inq_varid(ncid, sdsname, &sds_id);
		
            /* Read the data. */
            if (nc_get_var(ncid, sds_id, &ssttmp[0][0]) !=0){
                fHY1C_out(stderr,"-E- %s line %d:  Error reading %s from %s.\n",
                __FILE__,__LINE__,sdsname,sstfile);
		exit(1);
            }

            if (nc_get_att_float(ncid, sds_id, "scale_factor", &slope) !=0){
                fHY1C_out(stderr,"-E- %s line %d: error reading scale factor.\n",
		__FILE__,__LINE__);
		exit(1);
            }

            if (nc_get_att_float(ncid, sds_id, "add_offset", &offset) !=0){
	        fHY1C_out(stderr,"-E- %s line %d: error reading scale offset.\n",
                __FILE__,__LINE__);
                exit(1);
             }
            /* Close the file */
            if (nc_close(ncid) != 0){
                fHY1C_out(stderr,"-E- %s line %d: error closing %s.\n",
                __FILE__,__LINE__,sstfile);
                exit(1);
             }

        } else {
	/* Give HDF a whirl if netCDF fails... */

            sd_id = SDstart(sstfile, DFACC_RDONLY);
	    if (sd_id == FAIL){
		fHY1C_out(stderr,"-E- %s line %d: error reading %s.\n",
		__FILE__,__LINE__,sstfile);
		exit(1);
	    }

	    sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
	    status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
	    status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) &ssttmp[0][0]);
	    if (status != 0) {
		fHY1C_out(stderr,"-E- %s Line %d:  Error reading SDS %s from %s.\n",
		 __FILE__,__LINE__,sdsname,sstfile);
		exit(1);
	    } 
	    if (getHDFattr(sd_id,"scale_factor",sdsname,(VOIDP)&slope) != 0) {
		fHY1C_out(stderr,"-E- %s line %d: error reading scale factor.\n",
		__FILE__,__LINE__);
		exit(1);
	    }
	    if (getHDFattr(sd_id,"add_offset",sdsname,(VOIDP)&offset) != 0) {
		fHY1C_out(stderr,"-E- %s line %d: error reading scale offset.\n",
		__FILE__,__LINE__);
		exit(1);
            }
            status = SDendaccess(sds_id);
            status = SDend(sd_id);
	}

        /* rotate 180-deg and add wrapping border to simplify interpolation */
        /* new grid is -180.125,180.125 in i=0,361 and -90.125,90.125 in j=0,181    */

        for (j=0; j<ny; j++) {
 	    for (i=0; i<nx; i++) {
	        ii = (i < nx/2) ?  i+nx/2 : i-nx/2;
                if (ssttmp[j][i] > -999) 
                    sstref[j+1][ii+1] = ssttmp[j][i] * slope + offset;
                else
		    sstref[j+1][ii+1] = sstbad;
 	    }
            sstref[j+1][0]    = sstref[j+1][nx];
            sstref[j+1][nx+1] = sstref[j+1][1];
        }
        for (i=0; i<nx+2; i++) {
            sstref[0]   [i] = sstref[1][i];
            sstref[ny+1][i] = sstref[ny][i];
        }

        free(ssttmp);
    }


    /* locate LL position within reference grid */
    i = MAX(MIN((int) ((lon+180.0+dx/2)/dx),OI4NX+1),0);
    j = MAX(MIN((int) ((lat+ 90.0+dy/2)/dy),OI4NY+1),0);

    /* compute longitude and latitude of that grid element */
    xx = i*dx - 180.0 - dx/2;
    yy = j*dy -  90.0 - dy/2;

    /* bilinearly interpolate, replacing missing (land) values with average of valid values in box */
    t = (lon - xx)/dx;
    u = (lat - yy)/dy;

    ftmp = 0.0;
    ntmp = 0;
    if (sstref[j  ][i  ] > sstbad+1) {
	ftmp += sstref[j  ][i  ];
	ntmp++;
    }
    if (sstref[j  ][i+1] > sstbad+1) {
	ftmp += sstref[j  ][i+1];
	ntmp++;
    }
    if (sstref[j+1][i+1] > sstbad+1) {
	ftmp += sstref[j+1][i+1];
	ntmp++;
    }
    if (sstref[j+1][i  ] > sstbad+1) {
	ftmp += sstref[j+1][i  ];
	ntmp++;
    }
    if (ntmp > 0) {
	ftmp /= ntmp;
	reftmp[0][0] = (sstref[j  ][i  ] > sstbad+1 ? sstref[j  ][i  ]: ftmp);
	reftmp[0][1] = (sstref[j  ][i+1] > sstbad+1 ? sstref[j  ][i+1]: ftmp);
	reftmp[1][1] = (sstref[j+1][i+1] > sstbad+1 ? sstref[j+1][i+1]: ftmp);
	reftmp[1][0] = (sstref[j+1][i  ] > sstbad+1 ? sstref[j+1][i  ]: ftmp);

	sst = (1-t)*(1-u) * reftmp[0][0]
	    + t*(1-u)     * reftmp[0][1]
	    + t*u         * reftmp[1][1]
	    + (1-t)*u     * reftmp[1][0];
	
	/* apply skin correction */
	sst -= 0.17;

    } else
        sst = sstbad;

    return(sst);    
}



/* ----------------------------------------------------------------------------------- */
/* get_oisst() - read and interpolate flat binary Reynolds OIV2 SST                    */
/*                                                                                     */
/* The files were written in IEEE binary (big-endian). Each file contains four FORTRAN */
/* records described as follows:                                                       */
/*                                                                                     */
/*    rec 1: date and version number        (8 4-byte integer words)                   */
/*    rec 2: gridded sst values in degC     (360*180 4-byte real words)                */
/*    rec 3: normalized error variance      (360*180 4-byte real words)                */
/*    rec 4: gridded ice concentration      (360*180 1-byte integer words)             */
/*                                                                                     */
/* B. Franz, SAIC, May 2004.                                                           */
/* ----------------------------------------------------------------------------------- */

#define OINX 360
#define OINY 180

float get_oisst(char *sstfile, float lon, float lat)
{
    static int   firstCall = 1;
    static int   nx = OINX;
    static int   ny = OINY;
    static float dx = 360.0/OINX;
    static float dy = 180.0/OINY;
    static float sstref[OINY+2][OINX+2];

    float sst = BAD_FLT;
    int   i,j,ii;
    float xx,yy;
    float t,u;

    if (firstCall) {

        FILE *fp = NULL;
        float ssttmp[OINY][OINX];
        int32_t syear,smon,sday;
        int32_t eyear,emon,eday;
        int32_t ndays,version;

        firstCall = 0;

        if ((fp = fopen(sstfile,"r")) == NULL) {
            HY1C_out("Error opening SST reference file %s for reading.\n",sstfile);
            exit(1);
        }

        if (fseek(fp,4,SEEK_SET) < 0) {
            HY1C_out("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}          
        if (fread(&syear   ,sizeof(int32_t),1,fp) != 1) {
            HY1C_out("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}
        fread(&smon    ,sizeof(int32_t),1,fp);
        fread(&sday    ,sizeof(int32_t),1,fp);
        fread(&eyear   ,sizeof(int32_t),1,fp);
        fread(&emon    ,sizeof(int32_t),1,fp);
        fread(&eday    ,sizeof(int32_t),1,fp);
        fread(&ndays   ,sizeof(int32_t),1,fp);
        fread(&version ,sizeof(int32_t),1,fp);
        fseek(fp,4,SEEK_CUR);

        if ( endianess() == 1 ) {
            swapc_bytes((char *)&syear    ,sizeof(int32_t),1);
            swapc_bytes((char *)&smon     ,sizeof(int32_t),1);
            swapc_bytes((char *)&sday     ,sizeof(int32_t),1);
            swapc_bytes((char *)&eyear    ,sizeof(int32_t),1);
            swapc_bytes((char *)&emon     ,sizeof(int32_t),1);
            swapc_bytes((char *)&eday     ,sizeof(int32_t),1);
            swapc_bytes((char *)&ndays    ,sizeof(int32_t),1);
            swapc_bytes((char *)&version  ,sizeof(int32_t),1);
        }

        HY1C_out("Loading OI Reynolds SST reference from %s\n",sstfile);
        HY1C_out("  file start date (y/m/d): %d / %d / %d\n",syear,smon,sday);
        HY1C_out("  file end   date (y/m/d): %d / %d / %d\n",eyear,emon,eday);
        HY1C_out("  days composited: %d\n",ndays);
        HY1C_out("  file version: %d\n",version);
        HY1C_out("\n");

        if (fseek(fp,4,SEEK_CUR) < 0) {
            HY1C_out("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}          
        if (fread(&ssttmp[0][0],sizeof(float),nx*ny,fp) != nx*ny) {
            HY1C_out("Error reading SST reference file %s.\n",sstfile);
            exit(1);
	}
        fclose(fp);

        if ( endianess() == 1 )
            swapc_bytes((char *)&ssttmp[0][0],4,nx*ny);

        /* rotate 180-deg and add wrapping border to simplify interpolation */
        /* new grid is -180.5,180.5 in i=0,361 and -90.5,90.5 in j=0,181    */

        for (j=0; j<ny; j++) {
 	    for (i=0; i<nx; i++) {
	        ii = (i < nx/2) ?  i+nx/2 : i-nx/2;
                sstref[j+1][ii+1] = ssttmp[j][i];
 	    }
            sstref[j+1][0]    = sstref[j+1][nx];
            sstref[j+1][nx+1] = sstref[j+1][1];
        }
        for (i=0; i<nx+2; i++) {
            sstref[0]   [i] = sstref[1][i];
            sstref[ny+1][i] = sstref[ny][i];
        }
    }


    /* locate LL position within reference grid */
    i = MAX(MIN((int) ((lon+180.0+dx/2)/dx),OINX+1),0);
    j = MAX(MIN((int) ((lat+ 90.0+dy/2)/dy),OINY+1),0);

    /* compute longitude and latitude of that grid element */
    xx = i*dx - 180.0 - dx/2;
    yy = j*dy -  90.0 - dy/2;

    /* bilinearly interpolate */
    t = (lon - xx)/dx;
    u = (lat - yy)/dy;

    sst = (1-t)*(1-u) * sstref[j  ][i  ]
        + t*(1-u)     * sstref[j  ][i+1]
        + t*u         * sstref[j+1][i+1]
        + (1-t)*u     * sstref[j+1][i  ];

    /*
    sst = sstref[j][i];
    */

    return(sst);    
}

/* ----------------------------------------------------------------------------------- */
/* get_sstref() - retrieves reference sea surface temperature              .           */
/*                                                                                     */
/* B. Franz, SAIC, May 2004.                                                           */
/* ----------------------------------------------------------------------------------- */
#define PATHCLIM 1
#define OISSTBIN 2
#define AVHRRSST 3
#include "smi_climatology.h"
float get_sstref(char *sstfile, float lon, float lat, int day)
{
    float sst = sstbad;
    int32_t  sd_id;

    if (format < 0) {

        /* Does the file exist? */
        if (access(sstfile, F_OK) || access(sstfile, R_OK)) {
            HY1C_out("-E- %s: SST input file '%s' does not exist or cannot open.\n",
                   __FILE__, sstfile);
            exit(1);
        }

        /* What is it? */
        char title[255] = "";
        if (nc_open(sstfile, NC_NOWRITE, &sd_id) == 0) {
          /* Format is netCDF */
    	  if (nc_get_att_text(sd_id, NC_GLOBAL, "title", title) == 0);  	
	        if (strstr(title,"Daily-OI-V2") != NULL) {
                    format = AVHRRSST;
		} else {                                                        
                    HY1C_out("-E- %s : Unable to initialize SST file\n",__FILE__);
                    HY1C_out("%s %d\n",sstfile,day);                              
                    exit(1);                                                    
                }      
                nc_close(sd_id); 
        } else {
            sd_id = SDstart(sstfile, DFACC_RDONLY);
            if (sd_id != FAIL) {
	        /* Format is HDF */
	        if (SDreadattr(sd_id,SDfindattr(sd_id,"Title"),(VOIDP)title) == 0) {
	            if (strstr(title,"SST Climatology") != NULL) {
                        format = PATHCLIM;
                        if (smi_climatology_init(sstfile,day,SST) != 0) {
                            HY1C_out("-E- %s : Unable to initialize SST file\n",__FILE__);
                            HY1C_out("%s %d\n",sstfile,day);
                            exit(1);
		        }
		    }
	        }
                SDend(sd_id);             
            } else {

	        /* Format is not netCDF or HDF.  Assuming OISST Binary. */
  	        format = OISSTBIN;
            }
	}

        /* need a back-up climatology for NRT products */
        if (format != PATHCLIM) {
            char *tmp_str;
            char  file   [FILENAME_MAX] = "";
            if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
                HY1C_out("OCDATAROOT environment variable is not defined.\n");
               exit(1);
            }
            strcpy(file,tmp_str); strcat(file,"/common/sst_climatology.hdf"); 
            if (smi_climatology_init(file,day,SST) != 0) {
                HY1C_out("-E- %s : Unable to initialize SST file\n",__FILE__);
                HY1C_out("%s %d\n",file,day);
                exit(1);
            }
	}

    }

    switch (format) {
      case PATHCLIM:
        sst = smi_climatology(lon,lat,SST);
        break;
      case OISSTBIN:
        sst = get_oisst(sstfile,lon,lat);
        break;
      case AVHRRSST:
        sst = get_avhrrsst(sstfile,lon,lat);
        break;
      default:
        HY1C_out("-E- %s: unknown SST input file format for %s.\n",__FILE__, sstfile);
        break;
    }

    if (sst < (sstbad+1)) {
        sst = smi_climatology(lon,lat,SST);
    }

    return(sst);
}


